#ifndef PCPS_NUMBERS_HEADER
#define PCPS_NUMBERS_HEADER

#include <complex>

#include <src/exception.h>

namespace pcps {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   returns the number one in the appropriate type or raises an error if the type is not
  ///          supported
  ///
  /// \return the number one
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> inline S unity() {
    throw pcps::Exception( ( boost::format("no specialization of pcps::unity() for type \"%s\"")
                             % typeid(S).name() ).str() );
    return S();
  }
  // specialization for double
  template <> inline double unity<double>() { return 1.0; }
  // specialization for complex
  template <> inline std::complex<double> unity< std::complex<double> >() { return std::complex<double>(1.0, 0.0); }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   returns the number zero in the appropriate type or raises an error if the type is not
  ///          supported
  ///
  /// \return the number zero
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> inline S zero() {
    throw pcps::Exception( ( boost::format("no specialization of pcps::zero() for type \"%s\"")
                             % typeid(S).name() ).str() );
    return S();
  }
  // specialization for double
  template <> inline double zero<double>() { return 0.0; }
  // specialization for complex
  template <> inline std::complex<double> zero< std::complex<double> >() { return std::complex<double>(0.0, 0.0); }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   returns the square root of -1 in the appropriate type or raises an error if the type
  ///          is not supported
  ///
  /// \return the number i
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> inline S imaginary_unity() {
    throw pcps::Exception( ( boost::format("no specialization of pcps::imaginary_unity() for type \"%s\"")
                             % typeid(S).name() ).str() );
    return S();
  }
  // specialization for double
  template <> inline double imaginary_unity<double>() { return 0.0; }
  // specialization for complex
  template <> inline std::complex<double> imaginary_unity< std::complex<double> >() { return std::complex<double>(0.0, 1.0); }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   returns the complex conjugate of a number
  ///
  /// \param[in]     value    the number to be conjugated
  ///
  /// \return the complex conjugate of value
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> inline S conj(const S value) {
    return value;
  }
  // specialization for complex
  template <> inline std::complex<double> conj< std::complex<double> >(const std::complex<double> value) {
    return std::conj(value);
  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   returns the real part of a number
  ///
  /// \param[in]     value    the number whose real part we want
  ///
  /// \return the real part of value
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> inline double real(const S value) {
    return double(value);
  }
  // specialization for complex
  template <> inline double real< std::complex<double> >(const std::complex<double> value) {
    return std::real(value);
  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   returns the imaginary part of a number
  ///
  /// \param[in]     value    the number whose imaginary part we want
  ///
  /// \return the imaginary part of value
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> inline double imag(const S value) {
    return 0.0;
  }
  // specialization for complex
  template <> inline double imag< std::complex<double> >(const std::complex<double> value) {
    return std::imag(value);
  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   returns the square of the norm of a number
  ///
  /// \param[in]     x        the number whose square norm we want
  ///
  /// \return the square norm of x
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  inline double square_norm(const double x) {
    return x*x;
  }
  inline double square_norm(const std::complex<double> x) {
    return x.real() * x.real() + x.imag() * x.imag();
  }

}

#endif
